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We employ a novel algorithm using a quasi-exact embedded-cluster matching technique as mini- 
mization method within a genetic algorithm to reliably obtain numerically exact ground states of the 
Edwards- Anderson XY spin glass model with bimodal coupling distribution for square lattices of up 
to 28 X 28 spins. Contrary to previous conjectures, the ground state of each disorder replica is non- 
degenerate up to a global 0(2) rotation. The scaling of spin and chiral defect energies induced by 
applying several different sets of boundary conditions exhibits strong crossover effects. This suggests 
that previous calculations have yielded results far from the asymptotic regime. The novel algorithm 
and the aspect-ratio scaUng technique consistently give 6s = —0.308(30) and 9c = —0.114(16) for 
the spin and chiral stiffness exponents, respectively. 
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Since the suggestion of Edwards and Anderson (EA) 
to capture the essence of spin glass behavior in a class of 
simple lattice models thirty years ago [1], the quest for 
their understanding has spurred an enormous research 
effort [2]. EA considered the Hamiltonian 



n 



Si • Sj 



(1) 



with 0(n) spins Si on a regular lattice with quenched, 
random and frustrated nearest-neighbor interactions Jij . 
Although substantial progress has been made in recent 
years in understanding Ising and vector spin glasses in 
finite dimensions D, mostly by the development and ap- 
plication of sophisticated numerical techniques, we still 
lack an undisputed theory of the spin glass phase [2] . Due 
to its relative simplicity, by far the most work has been 
devoted to the Ising spin glass [2]. However, much less 
advance has been made on models with continuous spins 
which are often more relevant to real materials [2]. 

The properties of the spin glass phase in the EA 
model are described by a scaling theory of the associ- 
ated zero-temperature fixed point [3]. The correspond- 
ing renormalization-group (RG) picture considers the 
scaling of the width of the distribution of random cou- 
plings, PL{Jij), with the coarse-graining length scale L, 
J{L) ~ JL^% defining the spin stiffness exponent 9s- 
Depending on whether 9s > or 9s < 0, the spin glass 
phase is stable or unstable against thermal fluctuations, 
respectively. Following a suggestion by Banavar et al. 
and McMillan [4], the scaling of J{L) can be inferred 
from monitoring the dependence of the energy of droplet 
or domain- wall excitations induced by a change of bound- 
ary conditions (BCs), giving rise to the name "domain- 
wall RG" (DWRG) method. For cases where 9s < 0, and 
thus the spin-glass transition temperature Tg — 0, such 
as for the EA Ising model in two dimensions (2D) [2], 9s 
also determines the critical behavior with the spin glass 
correlation length diverging as ^ ~ j^-i^n fQj- T | 0, where 



Vs — ~^/9s [3]. Furthermore, unless exact ground-state 
degeneracies occur, as for the Ising model with bimodal 
P{Jij) [3], 9s is the only non-trivial exponent, while the 
critical exponent 77 is simply 2 — D when 9s < [3] . Con- 
sequently, 2D models offer a crucial test bench for our 
understanding of spin glasses at low temperatures. 

Twenty years of research since the original DWRG 
work of Morris et al. [5] have not been able to set- 
tle a number of persistent controversies concerning the 
ground-state properties of the 2D XY spin glass. Firstly, 
it has been suggested that the ground state may possess 
non-trivial extensive degeneracies when P{Jij) is a dis- 
crete bimodal distribution [6-8] . Secondly, it was realized 
early [9] that the rotational symmetry of the XY spin 
glass is accompanied by a Z2 symmetry originating from 
the difference between proper and improper 0(n) rota- 
tions [10]. It has been suggested that the resulting Ising- 
like chirality variables may decouple from the rotational 
degrees-of-freedom, leading to different critical behavior 
for the spin and chiral variables [11]. For D = 2, where 
Tg = 0, this would entail distinct spin and chiral stiffness 
exponents, 9s — —1/vs and 9c = — l/i'c, respectively. Fi- 
nally, and most noteworthy, previous Monte Carlo (MC) 
[6] and DWRG studies [5, 12, 13] have yielded rather in- 
consistent values for 9s- This might be partly explained 
by the difficulty in obtaining ground state configurations 
of the model. Parallel alignment of the spins to their local 
molecular fields hi — Jij Sj is a necessary condition 
for metastability of the system (1). However, due to the 
broad spectrum of an exponential number of metastable 
states, the resulting commonly used [5, 12] iterative spin 
quench algorithm [14] almost never yields a ground state 
configuration. Additionally, experience with the simpler 
2D Ising case shows that finite-size corrections as well as 
the dependence on the chosen pair of boundary condi- 
tions are generically large [15, 16]. Hence it seems likely 
that the observed inconsistencies for the XY model to 
date are due to system-size restrictions, improper finite- 
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FIG. 1: Histogram H{E) of minimum energies obtained from 
repeated runs using the spin-quench method (green, diago- 
nally hatched bars), simulated annealing (blue, solid bars) 
and the genetic embedded-matching technique (red, horizon- 
tally hatched bars) for a single disorder realization {Jij} of 
system-size 24 x 24. The insets show blow-ups of the region 
around the true ground state. The genetic matching was run 
here with a small population of A/ij ~ 64 replica; for Ao ~ 256, 
all runs converge to the rightmost bar of the top inset. 



size scaling analyses, and limitations in probing the true 
ground state behavior. 

Toulouse noted that the ground state problem for the 
2D Ising spin glass on a planar graph can be transformed 
to a minimization problem for the total length of energy 
strings on the dual lattice, connecting pairs of frustrated 
plaquettes, i.e., graph faces containing an odd number of 
negative bonds [17, 18]. It was later realized [19] that this 
constitutes a minimum- weight perfect matching problem, 
which is well known in graph theory and can be solved 
in polynomial time, such that the ground state of large 
2D Ising spin glass systems can be found exactly. In con- 
trast, the XY model ground state problem is seemingly 
not polynomial. However, in this paper we propose that 
a partial solution can be found in polynomial time by an 
embedding of Ising variables into the continuous spins, 
allowing us to obtain new results addressing the contro- 
versies alluded to above. This embedding is achieved by 
choosing a random direction r in spin space to decompose 
the spins a,s St ^ sf + ^ {S^ ■ r)r + S:^. A reflection 



.11 



Ri{r) of Si along the plane defined by r maps S, 
and S^ S:^. Hence, with respect to these local reflec- 
tions the Hamiltonian (1) decomposes asTi. = 

with ^'^■11 = -E(,,,)4<S''a"d 



Jij — Jij I 



r\\Sj-rl e[ = sign(S'j -r). 



(2) 



Thus, since the Ri{r) merely induce an inversion e[ — > 
— e[, the 0(n) model Hamiltonian (1) is formally iden- 
tified with that of an Ising model, if spin changes are 
restricted to the reflections Ri{r). One can then proceed 
as follows: decompose the 0(n) spins Si with respect 
to r and find the corresponding Ising ground state us- 
ing the matching technique [19]. This corresponds to a 



reflection of some of the Si and thus a new valid 0(n) 
model configuration. With Ti*"'"^ being invariant, this em- 
bedded ground state search decreases the total energy of 
(1) or leaves it constant. The full 0(n) symmetry can 
then be statistically recovered by sequential minimiza- 
tions for a series ri, r2, r3, . . . of random directions. We 
call this procedure "embedded matching". If (1) is in a 
ground state, all the embedded Ising systems must be in 
(one of) their respective ground state(s) as well. How- 
ever, stationarity of the process of successive embedded 
Ising-like matching minimization steps does not guaran- 
tee global minimum energy for the system (1) [20]. Thus, 
the corresponding artificial dynamics of exhibits metasta- 
bility, however with far less metastable states than the 
local spin-quench method [20]. For further improvement, 
and to find true ground states with high reliability, the 
embedded matching procedure is inserted as a minimiza- 
tion step into a specially tailored genetic algorithm [20]. 
Generally speaking, in a genetic algorithm, a popula- 
tion of A/q candidate ground-state configurations is being 
iteratively optimized by mixing or "crossing over" the 
"genetic material" of different candidate ground states 
and eliminating the less well adapted instances [21]. To 
achieve reasonable performance, this crossover opera- 
tion has to be chosen appropriately. Specifically, we 
are guided by the direct (visual) inspection of the spin 
configurations from different metastable states obtained 
by the embedded matching technique. There, due to 
the local spin rigidity, the predominant differences con- 
sist of (proper or improper) 0{n) rotations of rigid do- 
mains. Hence, to preserve the high level of optimiza- 
tion already obtained inside of domains at intermediate 
stages of the evolution, new offspring conflgurations are 
produced by randomly exchanging these (automatically 
determined) domains instead of single spins between the 
parent replica. Full details of the algorithm will be pre- 
sented elsewhere [20] (for a related method for the Ising 
case see Ref. [22]). This "genetic embedded- matching" 
(GEM) approach works very reliably already for small 
TVo as shown in Fig. 1. There, we compare the histograms 
of energies of metastable states found from statistically 
independent runs for the same ± J disorder conflguration 
{Jij} of a 24 X 24 system using either the simple spin- 
quench approach [14], the simulated annealing method, 
or the GEM technique. The first two methods give broad 
distributions of energies, whereas on this scale, the GEM 
always seems to yield the same energy, which is clearly 
below the range of energies regularly found by the other 
approaches. Only on examining the histograms at much 
higher resolution, do the GEM data get resolved into a 
small series of sharp peaks, corresponding to different en- 
ergy levels, cf. the upper inset of Fig. 1. On increasing A/o 
from A/q — 64, chosen for the runs in Fig. 1, to TVq — 256, 
the peaks displayed in the inset all collapse onto the peak 
of lowest energy on the right, corresponding to the true 
ground state. 
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FIG. 2: Log-log plot of the average domain- wall energies 
[|Ai5|]j for three sets of boundary conditions on square lat- 
tices as a function of system size L. The lines are fits of the 
form [|Ai5|]j ~ to the data. The black circles show the 
"random twist" result of Ref. [13] for comparison. Some data 
sets have been shifted vertically for better distinction. 



As a first result we find, perhaps surprisingly, that the 
ground states thus obtained for a given bimodal disor- 
der realization are not only identical in energy up to 
machine precision (15 digits) between statistically inde- 
pendent runs, but the final optimized spin configurations 
themselves are trivially related to each other by global 
0{n) transformations [20]. In other words, the ground 
states are unique and, in contrast to the bimodal Ising 
model [2], no accidental degeneracies occur. Hence, after 
averaging over disorder, the ground state is ordered and 
the spin correlation function is constant, implying 77 = 
[3] . This is in contrast to indications by MC simulations 
[6, 7] and Migdal-Kadanoff calculations [8], which pre- 
sumably did not probe the true ground state behavior. 

Strong dependence of the domain-wall scaling on the 
choice of BCs has been observed for the Ising spin glass 
[15, 23]. Further complications arise for the XY case 
due to the simultaneous presence of continuous (spin) 
and discrete (chiral) symmetries: for both periodic (P) 
and antiperiodic (AP) BCs, domain walls might be forced 
into the system due to the periodicity, such that the 
P / AP energy difference does not directly capture the en- 
ergy of single walls. In Ref. [13], it was attempted to 
alleviate this problem by introducing a twist along the 
boundary, which is included in the optimization process 
to yield "optimum twist" BCs. Nevertheless, additional 
chiral domain-walls might still occur in the measure- 
ments of spin domain- walls. In fact it has been found 
that for (quasi) one-dimensional XY systems both P/AP 
and reflective BCs asymptotically probe the chiral exci- 
tations [24]. These problems, resulting from a periodic 
constraint, can be circumvented by applying open and 
"domain- wall" (0/DW) BCs. Here, one ensures the in- 
sertion of single domain walls by comparing the ground 
state of a system with open BCs to one where the rela- 
tive orientations of spins linked across the boundary are 



either tilted by an angle tt for spin domain walls or re- 
flected along an arbitrary but flxed axis for chiral domain 
walls by the introduction of very strong bonds [20, 25]. 
In addition, and for comparison, we consider P/AP and 
random-antirandom (R/AR) BCs as well, the latter fix- 
ing the boundary spins in random relative orientations 
for one ground state computation (R) and in relatively 
TT-rotated orientations for AR BCs. In all cases, the 
edges with unaltered BCs are left open. Ground states 
were computed for systems of up to 28 x 28 spins, using 
5000 disorder realizations with = ±J at equal pro- 
portions. Figure 2 shows the results for the three sets 
of BCs together with fits of the asymptotically expected 
form [|Ai?|]j ^ to the data, where [•]j denotes the 
average over disorder. The results for P/AP BCs show 
a pronounced crossover from 6 ~ —0.724(21) for L < 12 
to 9 = —0.433(26) for L > 16, the first value being com- 
patible with the "random twist" data of Ref. [13] drawn 
for comparison {9^ = —0.76), which are representative 
of previous results for P / AP BCs and small system sizes 
[12]. On the other hand, 9 = -0.433(26) is closer to the 
"optimum twist" result of Ref. [13], designed to allevi- 
ate the problem of trapped domain walls. Note that the 
apparent crossover length is compatible with the length 
below which no metastability occurs and the system be- 
haves like a spherical spin glass [5, 26]. The other BCs 
yield less negative values already for smaller system sizes, 
resulting in 9s — —0.519(30) for the R/AR combination 
and 9s = -0.207(12) for the 0/DW BCs. The scaling of 
the chiral domain- wall energies from 0/DW boundaries 
yields an only slightly negative value 9c = —0.090(23). 

The above usage of multiple pairs of BCs reveals the 
presence of pronounced finite-size corrections, even for 
the already larger system sizes considered here compared 
to previous studies [5, 12, 13]. Part of these corrections 
are due to irrelevant scaling fields and sub-leading analyt- 
ical terms, giving rise to the general form [|Ai?|]j(L) — 
AL^ + BL-'^ + C/L + D/L"^ H . For a proper reso- 
lution of these contributions, much larger system sizes, 
out of the reach of current numerical methods, would be 
necessary. Thus, we have to restrict ourselves here to a 
successive omission of data points from the small-L side 
to extrapolate towards L — > 00. Additional corrections, 
however, result from the dependence on the considered 
pair of BCs. For the Ising system, it has been argued 
that such corrections might be suppressed by consider- 
ing L X M systems (the change of BCs happening along 
the edges of length L) with aspect ratios R = M/L ^ 1 
[23]. Neglecting for the time being the corrections listed 
above, the asymptotic scaling of defect energies should 
then follow the form [\AE\]j{L, M) = L'^F{R) with some 
scaling function F. In general, F{R) depends on the BCs 
applied [23]. However, there is no dependence on BCs 
for one-dimensional systems [3, 24], such that F{R) is 
independent of BCs in the limit R 00 and the corre- 
sponding corrections should disappear as more and more 
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FIG. 3: Aspect-ratio scaling of the stiffness exponents 9s and 
Oc for aspect ratios R = 1, 2 and 6 as a function of 1/R. 
The bottom data set corresponds to fits for the fixed R data 
restricted to L < 10 (see text). The solid lines show fits of 
the functional form 0{R) = 9{R = oo) + An/ R to the data. 



to be 9c = —0.114(16), clearly different from the spin ex- 
ponent , indicating spin-chirality decoupling, and close 
to the value 0^ = found for the bimodal Ising spin glass 
[25]. Yet, no exact degeneracies as occurring in the lat- 
ter case are found here. It would be very interesting to 
see whether the XY spin glass with Gaussian couplings 
shows a different behavior than the ±J case considered 
here. The 2D Heisenberg spin glass is another exciting 
problem to explore. 
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elongated systems are being considered. To investigate 
this, we determined the ground states of 5000 disorder 
replica and L = 4, 6, . . . , 16 for i? = 2 and L = 3, 4, . . . , 9 
for i? = 6 in addition to the data for the square sys- 
tems [R = \) for the different sets of BCs. Figure 3 
shows the estimated stiffness exponents as a function of 
R for P/AP and 0/DW BCs together with fits of the 
functional form 9{R) = 0{R = oo) + Ar/R to the data, 
which is inspired by the results for the Ising case [15, 23]. 
The scaling corrections at fixed R listed above are taken 
into account by including only the largest lattice sizes 
in the fits of [|Ai?|],7 ^ . For comparison, the bot- 
tom dataset of Fig. 3 shows the results from including 
a fixed range of sizes L < 10 for each aspect ratio R, 
thus admixing the two correction effects. The fits re- 
sult in consistent asymptotic estimates of the spin stiff- 
ness exponent of 6'/(i? = oo) = -0.338(20) from P/AP 
BCs and of 94R = oo) = -0.308(30) from 0/DW BCs, 
indicating that the asymptotic regime is indeed being 
probed. The chiral exponent dc, on the other hand, de- 
pends only weakly on R, and the asymptotic estimate 
dc{R = oo) = —0.114(16) is clearly different from 9s- 

In conclusion, we have developed a novel quasi-exact 
algorithm to determine the ground state of 2D 0(n) spin 
glasses. Considering for specificity the 2D XY spin glass 
model with bimodal distribution of random exchange 
couplings Jij , we have shown from computations for rel- 
atively large systems sizes that, as argued in Ref. [13], 
defect- wall calculations from P/AP BCs indeed suffer 
from large finite-size corrections due to the periodic con- 
straint. Using aspect-ratio scaling, however, they are 
found to asymptotically yield the same scaling behav- 
ior as the less ambiguous 0/DW BCs showing less pro- 
nounced corrections, and we quote the latter result as our 
final estimate, 6s = —0.308(30). This might be compared 
with 0s ~ —0.28 for the 2D Ising case with Gaussian cou- 
pling distribution [15, 25]. The chiral exponent is found 
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